function [ top1Ts, top3Ts, objTr, objTs,time, k] = R1MP_f( X_tr, Y_tr, X_ts, Y_ts, maxRank )

%R1MP Rank-one matrix pursuit full observation 
% Input: Observations of the label matrix y and its properties, sample matrix X 
% Output: Low-rank approximation of the feature matrix W


tol = 1e-3;
X_tr = X_tr';
% X_tr = sparse(X_tr);
[N,L] = size(Y_tr);
D = size(X_tr,1);
Y_trorig = Y_tr;
%%%%%%%%%%%%%%add missing%%%%%%%%%%%%%%%%%%
% missRate = 0.8;
% Y_tr(find((rand(N,L)<missRate))) = -1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[Ik,Jk] = find(Y_tr~=-1);
known = find(Y_tr~=-1);
y = Y_tr(known);
size_known = numel(known);
% [Ik, Jk] = ind2sub([N,L],known);
objTr = zeros(maxRank,1);
objTs  =zeros(maxRank,1);
top1Ts = zeros(maxRank,1);
top3Ts = zeros(maxRank,1);
time = zeros(maxRank,1);
r = ones(maxRank,1);
U = [];
V = [];
temp = [];
% XtM = [];

yy = [];
 x = randn(D,1);
% x = randn(L,1);
k = 0;
XtW = zeros(size_known,1);

while (k<maxRank && r(k+1)>tol)
    tstart = tic;
    %Find the lesft and right singular value of current residual
    R = y - XtW;
    r(k+1) = sqrt(sum(R.^2)/size_known);
    R = sparse(Ik,Jk,R,N,L); %%%%%when has missing
%      R = reshape(R,N,L);  %%%dealing with full labels
%     if k >0
%         x = x - U*(U'*x);
%     end
    [u1, v1] = powerM(X_tr,R,x,100,tol); 

%     m1 = u1*v1';
%     m1 = m1(:);

    
    temp1 = X_tr'*u1;
    XtM1 = temp1(Ik).*v1(Jk);  %XtM1 = vec_o(Xt*Mk)=vec(Xt*u1*v1')

    % update coefficients
    yy = [yy;XtM1'*y];
    
    
    if k ~= 0
        bb = zeros(k,1);
        for  i = 1:k
            utemp = temp(Ik,i);
            vtemp = V(Jk,i);
            bb(i) = (utemp.*vtemp)'*XtM1;
        end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%any better solution?%%%%%%%%%%%

        cc = XtM1'*XtM1;
        Minv = inverse_incremental( Minv, bb, cc );
    else
        Minv = 1/(XtM1'*XtM1);
    end
    c = Minv*yy; 
    U = [U u1];
    V = [V v1];
    temp = [temp temp1];
%     XtM = [XtM XtM1];
   
%   m = [m m1];

    % update current learned matrix
    XtW = temp*diag(c)*V';
    XtW = XtW(known);
    k = k+1;
    temptime(k) = toc(tstart);
    
    W_est = U(:,1:k)*diag(c)*V(:,1:k)';
    %%%Evaluation on the training data%%%%%%%%%%%
    %%TOP 1 RATE AND TOP3 RATE
    Y_estTr = X_tr'*W_est;
    objTr(k) = sum((Y_tr(known)-Y_estTr(known)).^2);
    % [top1Tr, top3Tr] = topRate(Y_estTr, Y_trorig);
    % % aucTr = AUC(Y_estTr, Y_tr);

    %%%Evaluation on the test data%%%%%%%%%%%
    %%TOP 1 RATE AND TOP3 RATE
    Y_estTs = X_ts*W_est;
    [top1Ts(k), top3Ts(k)] = topRate(Y_estTs, Y_ts);
    objTs(k) = sum((Y_ts(:)-Y_estTs(:)).^2);
    %  aucTs = AUC(Y_estTs, Y_ts);
     if k>1
        if abs(objTr(k)-objTr(k-1))<1e-3
            break;
        end
     end
    
     if k>3
        if objTs(k)>objTs(k-1) && objTs(k-1)>objTs(k-2) && objTs(k-2)>objTs(k-3)
            break;
        end
    end
    
    
end
for i = 1:k
    time(i) = sum(temptime(1:i));
end


plot([5:1:k],log(r(5:k)));
xlabel('number of iterations');
ylabel('RMSE in logarithm scale');
title('Bibtex - full observation');
print -dtiff -r500 mr0

end




